IMAGE DECONVOLUTION UNDER POISSON NOISE USING SPARSE 
REPRESENTATIONS AND PROXIMAL THRESHOLDING ITERATION 

F.-X. Dupe", M.J. Fadili" and J.-L. Starck^ 

" GREYC UMR CNRS 6072 DAPNIA/SEDI-SAP CEA-Saclay 
14050 Caen France 91191 Gif-sur-Yvette France 



ABSTRACT 

We propose an image deconvolution algorithm wiien the data is 
contaminated by Poisson noise. The image to restore is assumed 
to be sparsely represented in a dictionary of waveforms such as 
the wavelet or curvelet transform. Our key innovations are: First, 
we handle the Poisson noise properly by using the Anscombe 
variance stabilizing transform leading to a non-linear degrada- 
tion equation with additive Gaussian noise. Second, the decon- 
volution problem is formulated as the minimization of a convex 
functional with a data-fidelity term reflecting the noise proper- 
ties, and a non-smooth sparsity-promoting penalties over the im- 
age representation coefficients (e.g. fi-norm). Third, a fast iter- 
ative backward-forward splitting algorithm is proposed to solve 
the minimization problem. We derive existence and uniqueness 
conditions of the solution, and establish convergence of the it- 
erative algorithm. Experimental results are carried out to show 
the striking benefits gained from taking into account the Pois- 
son statistics of the noise. These results also suggest that using 
sparse-domain regularization may be tractable in many deconvo- 
lution applications, e.g. astronomy or microscopy. 

Index Terms — Deconvolution, Poisson noise. Proximal 
iteration, forward-backward splitting. Iterative thresholding. 
Sparse representations. 

1. INTRODUCTION 

Deconvolution is a longstanding problem in many areas of sig- 
nal and image processing (e.g. biomedical imaging 1 1 1, astron- 
omy |2|, remote-sensing, to quote a few). For example, research 
in astronomical image deconvolution has recently seen consider- 
able work, partly triggered by the Hubble space telescope (HST) 
optical aberration problem at the beginning of its mission. In 
biomedical imaging, researchers are also increasingly relying on 
deconvolution to improve the quality of images acquired by con- 
focal microscopes. Deconvolution may then prove crucial for 
exploiting images and extracting scientific content. 

There is an extensive literature on deconvolution problems. 
One might refer to well-known dedicated monographs on the 
subject. In presence of Poisson noise, several deconvolution 
methods have been proposed such as Tikhonov-Miller inverse 
filter and Richardson-Lucy (RL) algorithms; see 1 1 2 1 for an ex- 
cellent review. The RL has been used extensively in applications 
because it is adapted to Poisson noise. The RL algorithm, how- 
ever, amplifies noise after a few iterations, which can be avoided 
by introducing regularization. In |3|, the authors presented a 
Total Variation (TV)-regularized RL algorithm, and Starck et al. 
advocated a wavelet-regularized RL algorithm 0. 



In the context of deconvolution with gaussian white noise, 
sparsity-promoting regularization over orthogonal wavelet co- 
efficients has been recently proposed (4l|5 |. Generalization to 
frames was proposed in (6l|7|. In (8), the authors presented an 
image deconvolution algorithm by iterative thresholding in an 
overcomplete dictionary of transforms. However, all sparsity- 
based approaches published so far have mainly focused on 
Gaussian noise. 

In this paper, we propose an image deconvolution algorithm 
for data blurred and contaminated by Poisson noise. The Poisson 
noise is handled properly by using the Anscombe VST, leading to 
a non-linear degradation equation with additive Gaussian noise, 
see (|2}. To regularize the solution, we impose a sparsity prior 
on the representation coefficients of the image to restore, e.g. 
wavelet or curvelet coefficients. Then, the deconvolution prob- 
lem is formulated as the minimization of a convex functional with 
a data-fidelity term reflecting the noise properties, and a non- 
smooth sparsity-promoting penalties over the image representa- 
tion coefficients. Inspired by the work in ^5], a fast proximal 
iterative algorithm is proposed to solve the minimization prob- 
lem. We also provide an analysis of the optimization problem 
and establish convergence of the iterative algorithm. Experimen- 
tal results are carried out to compare our approach and show the 
striking benefits gained from taking into account the Poisson na- 
ture of the noise. 

Notation 

Let Ti a real Hilbert space, here a finite dimensional vector sub- 
space of R" . We denote by 1 1 . 1 1 j the norm associated with the in- 
ner product in 7i, and I is the identity operator on Ti. x and a are 
respectively reordered vectors of image samples and transform 
coefficients. A function / is proper if it is not identically +oo 
everywhere. A function / is coercive, if lim||^||^^+oo ./ [x) = 
+00. ro(H) is the class of all proper lower semi-continuous 
convex functions from H to ] — cxj, +00]. The subdifferential of 
/ is denoted df. 

2. PROBLEM STATEMENT 

Consider the image formation model where an input image x is 
blurred by a point spread function (PSF) h and contaminated by 
Poisson noise. The observed image is then a discrete collection 
of counts y = {yi)i^i^n where n is the number of samples. 
Each count yi is a realization of an independent Poisson random 
variable with a mean {h®x)i, where (D is the circular convolution 
operator. Formally, this writes yi ^ V {{JuS) x)i). 

A naive solution to this deconvolution problem would be to 
apply traditional approaches designed for Gaussian noise. But 



this would be awkward as (i) the noise tends to Gaussian only 
for large mean (/i ® x)^ (central limit theorem), and (ii) the noise 
variance depends on the mean anyway. A more adapted way 
would be to adopt a bayesian framework with an adapted anti- 
log-likelihood score reflecting the Poisson statistics of the noise. 
Unfortunately, doing so, we would end-up with a functional 
which does not satisfy some key properties (the Lipschitzian 
property in |5|), hence preventing us from using the backward- 
forward splitting proximal algorithm to solve the optimization 
problem. To circumvent this difficulty, we propose to handle the 
noise statistical properties by using the Anscombe VST defined 
as 

z^ = 2y^^r+|^, l^i^n. (1) 

Some previous authors |9| have already suggested to use the 
Anscombe VST, and then deconvolve with wavelet-domain reg- 
ularization as if the stabilized observation Zi were linearly de- 
graded by h and contaminated by additive Gaussian noise. But 
this is not valid as standard asymptotic results of the Anscombe 
VST state that 

= 2^(h®x)^ + l+e, e~AA(0, 1) (2) 

where e is an additive white Gaussian noise of unit variance. In 
words, z is non-linearly related to x. In Section 1431 we pro- 
vide an elegant optimization problem and a fixed point algorithm 
taking into account such a non-linearity. 

3. SPARSE IMAGE REPRESENTATION 

Let X £ Ti. be m ^/n x -y/n image, x can be written as the 
superposition of elementary atoms ip^ parameterized by 7 G T 
such that: 

X = a-yip-f = "1>Q, \T\ = L (3) 
7ex 

We denote by <I> the dictionary i.e. the n x L matrix whose 
columns are the generating waveforms {'fij)^^^^ all normalized 
to a unit i'2-norm. The forward transform is then defined by a 
non-necessarily square matrix T = G R^^" with L ^ n. 
When L > n the dictionary is said to be redundant or over- 
complete. In the case of the simple orthogonal basis, the inverse 
transform is trivially <I> = T^. Whereas assuming that T is a 
tight frame implies that the frame operator satisfies T^T = AI, 
where A > is the tight frame constant. For tight frames, the 
pseudo-inverse reconstruction operator reduces to A~^T. 

Our prior is that we are seeking for a good representation 
of X with only few significant coefficients. This makes sense 
since most practical signals or images are compressible in some 
transform domain (e.g. wavelets, curvelets, DCT, etc). These 
transforms generally correspond to an orthogonal basis or a tight 
frame. In the rest of the paper, "I> will be an orthobasis or a tigth 
frame of Ti. 

4. SPARSE ITERATIVE DECONVOLUTION 

We first define the notion of a proximity operator, which was 
introduced in 1101 as a generalization of the notion of a convex 
projection operator. 

Definition 1 (MoreaufTOl). Let(p £ FoiH). Then, for every x G 
Ti., the function y ^ </'(j/) + 11^ ~ V\\^ /2 achieves its infimum 
at a unique point denoted by prox^ x. The operator prox^ : 
Ti. ^ Ti. thus defined is the proximity operator of ip. We define 
the reflection operator rprox^ = 2 prox^ —I. 



4.1. Optimization problem 

The class of minimization problems we are interested in can be 
stated in the general form |5]: 

argmin/i(a::) + /2(a;). (4) 

where /i G ro(H), /2 G ro(H) and /i is differentiable with 
K-Lipschitz gradient. We denote by M the set of solutions. 
From ^ and JS}, we immediately deduce the data fidelity term 

FoH o<i> (a), with (5) 
F:V^J2fM, /(r;0 = i(^.-2\/r?, + |) , 

where H denotes the (block-Toeplitz) convolution operator. 
From a statistical perspective, ^ corresponds to the anti-log- 
likelihood score. 

Adopting a bayesian framework and using a standard max- 
imum a posteriori (MAP) rule, our goal is to minimize the fol- 
lowing functional with respect to the representation coefficients 

Q 

(Pa.v-) : arg min J(q!) (6) 

a 

L 

J : a t-* F o H o ^ (a) + tc o ^ (a) + X ^ i^{cti), 
■' 1=1 

V ' 

/2(q) 

where we implicitly assumed that (ai)i^i^L are independent 
and identically distributed with a Gibbsian density oc e~'^^^°''K 
Notice that /2 is not smooth. The penalty function is chosen to 
enforce sparsity, A > is a regularization parameter and ic is the 
indicator function of a convex set C. In our case, C is the positive 
orthant. We remind that the positivity constraint is because we 
are fitting Poisson intensities, which are positive by nature. We 
have the following. 

Proposition 1. 

• /i is convex function. It is strictly convex if ^ is an ortho- 
basis and ker [H) = ^ (i.e. the spectrum of the PSF does 
not vanish ). 

• /i is continuously differentiable with a K-Lipschitz gradi- 
ent where 

K^{lfUA\\H\\l\\z\\^<+<^. (7) 

• (Px.ip) is a particular case of problem 

4.1.1. Characterization of the solution 

Proposition 2. Since J is coercive and convex, the following 
holds: 

1. Existence: (P\.^) has at least one solution, i.e. M 7^ 0. 

2. Uniqueness: (Pa,</)) has a unique solution if ^ is a basis 
and ker [H) = 0, or ifip is strictly convex. 

Proof: The existence is obvious because J is coercive. If $ 
is an ortho-basis and ker [H) = then /i is strictly convex and 
so is J leading to a strict minimum. Similarly, if is strictly 
convex, so is /2, hence J. 



4.1.2. Proximal iteration 

For notational simplicity, we denote by *!/ the function a 
4>{ai). The following useful lemmas are first stated: 

Lemma 1. The gradient of \/ fi is 

V/i(a) = o //* o VF o o $ (a) 

with 



VF(77) 



+ 2 



(8) 



(9) 



The proof is straightforward. 

Lemma 2. Let $ an orthobasis or a tight frame with constant c. 

1. If a G C then prox^ (a) — prox_;^^ (a), C' — {a\<^a G 
C}. 

2. Otherwise, let J]]^ //t(l — z/t) = +oo. Tafe 7*^ G 7i, anrf 



define the sequence of iterates: 
ft I rprox 



t+i t 

7 =7 



rprox,^, -I^ (7'), 



(10) 



where 
prox 



Vc = Pi-ox,^, = c"^<I>'^ o Pc o $ + (I - c ^^"^ o $) 
anrf "Pc w projector onto the positive orthant [Vc v)i ~ 
max(?7i, 0). Then, 



7' J and prox^^(a) = Vci-j)- 



(11) 



We are now ready to state our main proximal iterative algorithm 
to solve the minimization problem (Px,^): 

Theorem 1. For t ^ 0, let (pt)t be a sequence in ]0, +cx3[ such 
thatO < inftpt s; sup^ /it < {^^^^ / {2A\\H\\l\\z 
ao G Ti., for every t ^ 0, set 



Fix 



it+i = prox^t/2 (at - pt (V/i(Qt))) 



(12) 



where V/i and prox 



given by Lemma Q] and |2] Tfeen 
(ot)t^o converges (weakly) to a solution o/(Pa,i/))- 

Proof: We give a sketch of the proof. The main theorem on 
the proximal iteration can be found in [5, Theorem 3.4]. Hence, 
combining this theorem with Lemma [T] Lemma |2] and Proposi- 
tion[T] the result follows. 

■ 

Note that if the PSF h is low-pass normalized to a unit sum, then 
\\H\\l = ^- 

We now turn to prox^^ which is given by the following result: 

Theorem 2. Suppose that ip satisfies, (i) ip is convex even- 
symmetric , non-negative and non-decreasing on [0, +00), and 
'i/>(0) = 0. (ii) 4' w twice differentiable on R \ {0}. (Hi) ^ is 
continuous on R, it is not necessarily/ smooth at zero and admits 
a positive right derivative at zero tp^ (0) = lim^^g+ -^^^^ > 0. 
Then, the proximity operator pT0Xg^{/3) — a(/3) has exactly 
one continuous solution decoupled in each coordinate f3i : 



A proof of this theorem can be found in |I12I|. Among the 
most popular penalty functions tp satisfying the above require- 
ments, we have "(/'(q^O ~ \cti\, in which case the associated prox- 
imity operator is soft-thresholding. In this case, iteration l ll2t is 
essentially an iterative thresholding algorithm with a positivity 
constraint. 



5. EXPERIMENTAL RESULTS 

The performance of our approach has been assessed on several 
2D datasets, from which we here illustrate two examples. Our 
algorithm was compared to RL without regularization, RL with 
multi-resolution support wavelet regularization [2, RL-MRS], 
the naive proximal method that would assume the noise to be 
Gaussian (NaiveGauss), and the approach of |9| (AnsGauss). For 
all results presented, each algorithm was run with 200 iterations, 
except RL which was stopped when its MSE was smallest. 

In Fig[T] the original Lena image with a maximum intensity 
of 30 is depicted in (a), its blurred and blurred-l-noisy versions 
are in (b) and (c). With Lena, and for NaiveGauss, AnsGauss 
and our approach, the dictionary $ contained the curvelet tight 
frame. The deconvolution results are shown in Fig[nd)-(h). As 
expected, the RL is the worst as it lacks regularization. There are 
also noticeable artifacts in NaiveGauss, AnsGauss and RL-MRS. 
Our deconvolved image appears much cleaner. This visual im- 
pression is confirmed by quantitative measures of the quality of 
deconvolution, where we used both the mean ^i-error (adapted to 
Poisson noise), and the well-known MSE criteria. The mean £1- 
error for Lena is shown in Tab[T](similar results were obtained for 
the MSE). In general, our approach performs very well. At low 
intensity levels, RL-MRS has the smallest error very compara- 
ble to our approach. For the other intensity levels, our algorithm 
provides the best performance. At high intensity levels, Naive- 
Gauss is competitive. This comes as no surprise since this is an 
intensity regime where Poisson noise approaches the Gaussian 
behavior. On the other hand, the results reveal that AnsGauss 
performs poorly just after RL, probably because it does not han- 
dle properly the non-linearity of the degradation model (O after 
the VST. 

We further illustrate the capabilities of our approach on a 
confocal microscopy simulation. We have created a phantom of 
an image of a neuron dendrite segment with a mushroom-shaped 
spines, see Fig[2] The experimental settings were the same as 
for Lena except that the dictionary here contained the wavelet 
transform. The findings are similar to those of Lena both visually 
and quantitatively. 





Intensity regime 


Method 


^ 5 


^ 30 


100 


255 


Our method 


0.39 


0.93 


2.63 


7.21 


NaiveGauss 


0.59 


1.65 


3.56 


6.9 


AnsGauss 


0.87 


2.33 


4.61 


8.35 


RL-MRS 


0.35 


1.76 


4.31 


9.5 


RL 


1.97 


5.07 


9.53 


15.68 



Table 1. Mean ^i-error of all algorithms as a function of 
the intensity level. 




(d) (e) (f) 




(g) (h) 

Fig. 1. Deconvolution of Lena, (a) Original, (b) Blurred, 
(c) Blurred+noisy, (d) RL, (e) NaiveGauss, (f) AnsGauss, 
(g) RL-MRS, (h) Our algorithm. 

6. CONCLUSION 

In this paper, we presented a sparsity-based fast iterative tliresli- 
olding deconvolution algorithm that takes account of the pres- 
ence of Poisson noise. A careful theoretical characterization 
of the algorithm and its solution is provided. The encourag- 
ing experimental results clearly confirm the capabilities of our 
approach. 
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